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ABSTRACT 


In  spite  of  extensive  prior  research  on  the  generation  of  seismic  waves  by  underground  nuclear  explosions,  it  is  still 
not  possible  to  provide  a  complete  explanation  for  the  observed  wavefields,  particularly  at  regional  distances. 
Spherically  symmetric  explosion  models  embedded  in  layered  elastic  media  effectively  model  the  P  phases 
generated  by  explosions  and  the  major  characteristics  of  some  reflected  and  transmitted  phases.  Nonlinear 
axisymmetric  finite  difference  calculations  of  explosions,  including  gravity  and  the  effect  of  the  free  surface,  can 
model  a  more  realistic  explosion  source  that  directly  generated  shear  waves.  These  models  explain  more 
characteristics  of  explosion-generated  seismic  waves,  including  some  aspects  of  regional  shear  phases.  However,  it 
is  clear  that  linear  and  nonlinear  near-source  3D  effects  are  important  in  many  cases.  SH  waves  are  commonly 
observed  within  a  few  kilometers  of  explosions,  too  close  to  have  been  generated  by  (simple)  conversion  of  vertical 
and  radial  components,  and  often  larger  than  those  components.  Furthermore,  it  has  not  been  established  what 
impact  3D  effects  have  on  discriminants  and  on  explosion  yield  estimates.  It  is  important,  therefore,  to  be  able  to 
model  and  understand  how  3D  source  and  source  region  heterogeneity  affect  the  seismic  wavefield  and  what  impact 
this  has  on  parameters  used  for  nuclear  monitoring. 

In  this  new  project,  we  are  implementing  a  technique  that  allows  us  to  propagate  the  results  of  near-source  3D  finite 
difference  calculations  to  regional  and  teleseismic  distances,  and  to  use  the  results  to  investigate  the  impact  of  3D 
near-source  effects  on  regional  and  teleseismic  phases,  focusing  in  particular  on  the  generation  of  SH  phases  by 
explosion  sources.  Our  approach  is  to  perform  3D  explosion  source  region  calculations,  and  then  to  propagate  the 
wavefield  to  local,  regional  and  teleseismic  distances  using  layered  earth  Green’s  functions.  We  are  interested  in 
near-source  heterogeneities  in  both  the  nonlinear  and  linear  regimes,  and  therefore  require  both  nonlinear  and  linear 
3D  codes  to  model  the  source  region.  In  previous  projects,  we  have  used  two  nonlinear  codes,  STELLAR  and 
CRAM,  and  the  3D  linear  elastic  code  TRES3D.  STELLAR  is  an  Eulerian  2D  and  3D  code  well  suited  to 
performing  near-source  calculations,  and  CRAM  is  a  2D  axisymmetric  Lagrangian  code  that  is  capable  of  modeling 
the  transition  zone  between  the  vaporized  explosion  cavity  and  linear  elastic  region,  including  the  effects  of  gravity 
and  the  free  surface.  As  part  of  this  project,  we  are  developing  a  3D  Lagrangian  code  with  the  material  models  and 
capabilities  of  CRAM,  and  we  are  developing  the  interface  code  to  propagate  3D  near  source  waveforms  to  regional 
and  teleseismic  distances  using  full  waveform,  modal  and  far-field  body  wave  Green’s  functions. 

The  first  part  of  this  project  has  focused  on  the  development  of  these  numerical  tools  that  will  be  used  throughout 
the  project.  Our  plans  are  to  perform  calculations  for  realistic  near-source  3D  heterogeneity  in  the  nonlinear  region; 
to  propagate  these  calculations  to  local,  regional,  and  teleseismic  distances;  and  to  compare  results  with  explosion¬ 
generated  waveforms. 
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OBJECTIVES 

The  objective  of  this  project  is  to  investigate  the  generation  of  complex  seismic  waves  by  explosions  in  media  with 
3D  heterogeneity  using  a  method  based  on  the  exact  representation  theorem  for  propagating  complex  3D-source 
calculations  to  local,  regional,  and  teleseismic  distances  at  which  they  are  observed. 

RESEARCH  ACCOMPLISHED 

Introduction 


Source  physics,  near-source  scattering,  and  propagation  effects  are  all  important  to  understanding  seismic  phases 
used  in  nuclear  monitoring.  Significant  bodies  of  literature  exist  that  address  each  subject  individually.  In  addition  to 
extending  our  understanding  of  source  physics  and  near-source  scattering  to  include  the  effects  of  realistic  3D 
heterogeneity,  this  new  project  links  the  progression  of  energy  from  its  generation  by  the  source,  through  the 
near-source  region,  and  on  into  its  partitioning  among  local,  regional,  and  teleseismic  phases.  Distinguishing  the 
far-field  P  and  S  waves  enables  us  to  quantify  the  effects  of  3D  structure  on  both  P-  and  S-wave  generation. 
Complete  regional  waveforms  show  how  this  energy  is  partitioned  among  the  distinct  phases  important  to  event 
detection,  identification,  and  magnitude  estimation.  How  energy  is  distributed  among  surface  wave  modes 
determines  Lg  amplitudes  and  Rg  amplitudes,  including  near-source  conversions  between  these  phases.  Modal 
excitation  of  Lg  as  well  as  Rg  has  a  significant  depth  dependence  that  is  often  neglected  in  nuclear  monitoring 
studies  (Baker  et  al.,  2004),  and  these  amplitudes  can  be  disrupted  further  by  near  source  effects  and  source  region 
structure.  This  can  either  degrade  or  improve  discrimination  capability,  depending  on  how  well  it  is  understood. 

As  an  example  of  anomalous  (but  almost  universally  observed)  SH  waves  from  explosions,  we  present  data  from  the 
Soviet  Peaceful  Nuclear  Explosions  (PNE)  Quartz3  and  Quartz4  (Figure  1,  left).  Figure  1  (right)  shows  an  example 
of  strong  tangential  component  S  waves  from  the  Quartz  3  PNE.  This  8.5-kiloton  explosion  was  overburied  by  a 
factor  of  3  relative  to  normal  scaled  depth,  at  a  depth  of  759  meters.  This  event  should  therefore  be  a  good 
approximation  to  the  ideal  point,  spherical  explosion  source.  Nevertheless,  it  generated  strong  SH  waves.  Figure  2 
shows  an  image  of  the  source  location.  As  indicated  by  the  meandering  river,  the  site  is  quite  flat,  so  topography 
would  be  expected  to  have  little  effect  on  the  seismic  waves.  Tangential  component  S  at  large  distances  would  not 
be  surprising,  as  the  components  gradually  homogenize  with  distance;  however,  in  this  case  the  S  waves  exist 
strongly  across  the  entire  section,  suggesting  that  they  leave  the  source  region  as  SH  waves,  rather  than  being  rotated 
into  that  orientation  by  structure  during  propagation. 


Tangential  Sg  East  of  Quartz  3 


Figure  1.  Layout  of  stations  from  the  Quartz3  and  Quartz4  PNE  (left)  and  tangential  component  S  waves 

from  Quartz3  (right).  The  blue  part  of  the  time  series  and  the  red  lines  delimit  4  to  2  km/s.  Strong 
SH  waves  exist  over  the  entire  section. 
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Figure  2.  Quartz3  location  (left)  and  Quartz4  location  (right).  Each  image  is  24  km  by  20  km.  The  explosion 
location  is  in  the  center  of  each  figure.  Both  sites  are  quite  flat.  Granite  can  be  seen  at  the  surface  at 
the  Quartz  4  location. 


Quartz4  was  a  10-kt  explosion  at  slightly  shallower  depth  (557  m)  than  Quartz3  but  still  overburied  by  a  factor  of  2. 
The  explosion  is  in  granite,  and  granite  can  be  seen  at  the  surface  near  the  explosion  site  in  Figure  2.  Again,  there  is 
little  apparent  topography.  Kitov  (1997)  states  that  “Geological  structure  for  Quartz4  may  be  characterized  as  a 
homogeneous  half-space  covered  by  a  thin  (10  to  15  m)  layer  of  soil.”  Figure  3  shows  the  three  components  of  the 
Quartz4  seismograms  at  a  station  3 1  km  to  the  east.  There  is  a  strong  SH  wave  on  the  tangential  component  and  no 
obvious  trace  of  it  on  either  the  vertical  or  radial  component.  Very  similar  observations  were  made  for  the  deep  U.S. 
explosion  Rulison  (Figure  3,  right). 
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Figure  3.  Left:  Radial  (top),  transverse  (center),  and  vertical  (bottom)  seismograms  from  Quartz4  observed 
31  km  east  of  the  explosion.  Data  have  been  bandpass  Altered  between  1  and  5  Hz.  There  is  a 
strong  SH  phase  on  the  tangential  that  does  not  exist  on  the  radial  or  vertical  components. 

Right:  Similar  figure  showing  seismograms  at  a  distance  of  23  km  from  the  deep  U.S.  explosion 
Rulison  (from  Murphy  and  Archambeau,  1986). 

These  examples  (and  many  others)  point  to  the  need  to  understand  3D  source  effects.  While  it  is  not  surprising  to  see 
SH  waves  from  any  single  event,  what  is  surprising  is  that  they  are  always  present  and,  at  regional  distances,  scale 
with  yield  about  as  well  as  P  waves.  That  is,  while  there  are  many  mechanisms  such  as  near-source  scattering. 
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tectonic  release,  etc.,  that  can  generate  SH  waves,  all  of  these  effects  should  be  highly  variable  from  one  event  to  the 
next.  To  address  this  question,  we  need  to  look  at  the  types  of  3D  source  effects  that  can  exist,  the  range  of 
variability  that  would  be  expected  from  them,  and  whether  the  predictions  are  consistent  with  observations. 

Source  Region  Calculations  and  Propagation  Using  the  Representation  Theorem 

Our  approach  is  to  perform  3D  explosion  source  region  calculations  and  then  to  propagate  the  wavefield  to  local, 
regional,  and  teleseismic  distances,  using  layered-earth  Green’s  functions.  We  are  interested  in  near-source 
heterogeneities  in  both  the  nonlinear  and  linear  regimes  and  therefore  require  both  nonlinear  and  linear  3D  codes  to 
model  the  source  region.  In  previous  projects,  we  have  used  two  nonlinear  codes,  STELLAR  and  CRAM,  which  are 
described  briefly  in  Table  1  (we  have  also  used  the  ID  nonlinear  code  SKIPPER,  which  is  a  spherically  symmetric 
version  of  CRAM),  and  the  2D  and  3D  linear  elastic  code  TRES3D. 


Table  1.  Numerical  simulation  tools  to  be  used  in  this  project 


Numerical  Simulation  Tools 

STELLAR 

Eulerian  finite-difference  code.  Used  to  simulate  the  early  time 
history  of  the  explosion  shock.  It  handles  material  strength  correctly, 
which  is  difficult  for  an  Eulerian  code.  Uses  a  second-order  accurate 
Riemann  solver  scheme.  ID,  2D  planar  and  axisymmetric,  and  3D. 

CRAM 

Lagrangian  nonlinear  finite  difference  code.  Has  been  used 
extensively  for  axisymmetric  explosion  calculations.  A  3D  version 
of  the  code  is  being  developed  in  this  project. 

TRES3D 

Elastic  finite  difference  code.  2D  planar  and  axisymmetric,  and  3D. 

Elastodynamic  Representation  Theorem 

The  time-dependent  displacements  and  stresses  from  3D  source 
region  calculations  are  saved  on  a  monitoring  surface  located  outside 
the  region  of  nonlinear  response  and/or  3D  heterogeneity.  A 
numerical  implementation  of  the  representation  integral  is  then  used 
to  compute  the  corresponding  far-field  seismic  radiation. 

In  past  projects,  we  have  used  these  codes  in  the  following  ways: 


1 .  Axisymmetric  CRAM  was  used  together  with  the  representation  theorem  to  propagate  the  results  of 
nonlinear  axisymmetric  finite-difference  calculations  to  regional  and  teleseismic  distances  (Day  et  al., 

1987;  Stevens  et  al.,  1991;  Stevens  et  al.,  2004). 

2.  Axisymmetric  STELLAR  was  used  to  calculate  the  early  stages  of  an  explosion.  The  solution  at  an 
appropriate  time  was  overlaid  onto  CRAM  to  be  propagated  out  to  the  linear,  elastic  region.  The 
representation  theorem  was  used  to  propagate  the  waveform  to  regional  and  teleseismic  distances  (Rimer  et 
al.,  1994).  A  similar  technique  was  used  to  overlay  results  from  the  LANL  code  SOIL  onto  CRAM  and 
propagate  the  results  (Davis  et  al.,  1992). 

3.  3D  STELLAR  was  used  to  calculate  near-field  waveforms  from  explosions  in  rectangular  cavities  (Stevens 
et  al,  2006). 

4.  TRES3D  was  used  to  calculate  the  scattering  from  explosions  in  a  region  with  3D  heterogeneity  and 
topography  (Stevens  et  al.,  2004). 

In  the  current  project,  we  are  doing  the  following:  (1)  using  STELLAR  to  perform  very  near  source  nonlinear  3D 
calculations  and  TRES3D  to  perform  linear  elastic  near-source  and  source  region  3D  calculations,  (2)  developing  a 
3D  version  of  CRAM  to  model  explosions  from  the  source  out  through  the  very  important  nonlinear-to-linear 
transition  region,  and  (3)  completing  the  implementation  of  the  elastodynamic  representation  theorem  for  full 
waveforms,  modes,  and  body  waves. 

The  technique  for  propagating  numerical  calculations  using  the  representation  theorem  is  to  save  displacements  and 
stresses  on  a  monitoring  surface  surrounding  the  nonlinear  and/or  heterogeneous  region  of  the  calculation  and  then 
to  convolve  these  with  a  Green’s  function  appropriate  for  the  external  region  (Stevens  et  al.,  1991).  In  the  cases  that 
we  have  done  previously,  2D  axisymmetric  nonlinear  finite-difference  calculations  were  performed  to  model  the 
nuclear  explosion,  and  the  stresses  and  displacements  from  the  calculation  were  saved  on  a  cylindrical  surface  in  the 
elastic  region  outside  the  region  of  complex  nonlinear  behavior.  We  then  invoked  the  representation  theorem  and 
integrated  the  stresses  and  displacements  with  an  axisymmetric  Green’s  function  to  calculate  the  displacement  at  any 
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point  outside  the  calculation.  We  performed  such  calculations  in  2D,  using  Green’s  functions  for  far-field  body 
waves,  for  modes,  and  for  full  regional  waveforms  using  wavenumber  integration.  The  equations  for  the  Green’s 
functions  for  surface  waves  are  given  by  Bache  et  al.  (1982).  The  Green’s  functions  for  the  complete  seismograms 
are  computed  using  a  ring  load  source,  from  an  algorithm  based  on  the  work  of  Luco  and  Apsel  (1983)  and  Apsel 
and  Luco  (1983).  The  Green’s  functions  for  body  waves  are  generated  by  a  procedure  similar  to  that  described  by 
Bache  and  Harkrider  (1976),  using  a  saddle  point  approximation  to  calculate  a  far-field  plane  wave  for  a  given 
takeoff  angle  from  a  source  in  a  plane-layered  medium.  Our  objective  in  using  multiple  types  of  Green’s  functions  is 
to  gain  as  much  insight  as  possible  into  the  nature  of  the  seismic  wavefield  generated  by  the  source.  An  important 
part  of  the  current  project  is  to  adapt  these  techniques  to  propagate  seismic  waves  from  3D  source  calculations. 
Although  any  closed  surface  can  be  used  for  representation  theorem  integration,  we  use  a  cylindrical  surface  for 
axisymmetric  problems  and  a  rectangular  surface  for  3D  problems. 

3D  Implementation  of  the  Representation  Theorem 

The  key  to  extending  the  axisymmetric  representation  theorem  discussed  above  to  3D  is  to  recognize  that  while  the 
deformation  in  the  source  region  may  be  arbitrarily  complex,  if  the  structure  can  be  approximated  as  a  plane-layered 
medium  outside  of  the  source  region,  then  the  known  Green’s  function  for  a  plane-layered  medium  applies  (this  also 
assumes  that  we  can  neglect  the  interaction  of  any  backscattered  waves  returning  to  the  source  region  after  leaving 
it).  Note  that  the  representation  theorem  is  exact.  That  is,  no  matter  how  complex  the  3D  motion  is  on  the  source 
region  boundary,  it  will  be  correctly  propagated  by  the  representation  theorem.  The  following  benchmark  test 
demonstrates  the  performance  of  the  method  by  comparing  the  results  with  equivalent  finite  difference  calculations 
and  wavenumber  integration  seismograms. 

The  representation  theorem  states  that  displacement  at  an  observation  point  is  made  up  of  contributions  due  to  body 
forces  throughout  the  source  volume,  plus  contributions  due  to  the  traction  and  displacement  on  the  source  volume 
surface  (Aki  and  Richards,  1980).  In  the  three-dimensional  numerical  finite  difference  calculations,  we  save 
displacements  and  stresses  due  to  the  seismic  source  on  a  monitoring  surface  on  the  boundary  of  a  rectangle 
(5  planar  surfaces,  excluding  the  upper  surface),  and  calculate  Green’s  functions  from  each  point  on  the  monitoring 
surface  to  the  receiver.  Thus  the  synthetic  seismogram  at  the  receiver  point  X  outside  of  the  monitoring  surface  is 
obtained  by  integrating  over  the  monitoring  surface 

‘'‘■a/ 

in  the  frequency  domain,  where  G'  (^;  X)  and  5'^  (^;  X)  are  the  Green’s  function  and  the  stress  tensor  on  the 

monitoring  surface  due  to  a  unit  impulsive  force  at  X  in  direction  i,  is  the  traction  on  the  monitoring  surface 

due  to  the  seismic  source,  u  is  the  displacement  on  the  monitoring  surface,  and  n  is  the  normal  to  the  monitoring 
surface.  The  operator  *  denotes  convolution,  and  the  summation  convention  is  assumed. 

Equation  1  is  applicable  to  any  Green’s  function  for  the  exterior  model,  and  so  we  can  use  a  full  waveform  Green’s 
function,  far-field  body  wave  Green’s  fxinction,  and/or  modal  Green’s  function,  as  we  have  discussed  earlier  for 
axisymmetric  problems.  We  have  implemented  this  technique  for  full  waveform  seismograms,  calculating  the 
Green’s  functions  using  wavenumber  integration.  We  plan  to  implement  the  technique  for  body  waves  and  modes 
later  in  this  project. 

The  synthetic  seismograms  are  computed  using  the  following  steps; 

1 .  Displacements  and  stresses  are  saved  on  the  monitoring  surface  during  the  finite-difference  calculations. 

2.  If  necessary,  the  monitoring  solutions  are  resampled  onto  a  coarse  grid,  as  permitted  by  the  required 
resolution.  If  there  is  a  symmetry  boundary,  the  entire  monitoring  surface  is  constructed  first. 

3.  The  finite-difference  solutions  at  each  point  on  the  monitoring  surface  are  transformed  into  the  frequency 
domain. 

4.  The  displacement  and  stress  Green’s  functions  due  to  the  three  orthogonal  forces  at  the  receiver  location 
are  calculated  for  each  location  on  the  monitoring  surface  in  the  frequency  domain  (using  reciprocity). 

5.  Equation  ( 1 )  is  used  to  obtain  the  solution  at  the  receiver  in  the  frequency  domain. 
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6.  The  solution  is  transformed  back  to  the  time  domain. 

Figure  4  shows  an  example  of  a  test  of  this  technique.  In  this  test  case,  full  waveform  synthetics  are  generated  from 
a  3D  elastic  finite-difference  calculation  in  two  ways:  (1)  by  using  the  representation  theorem  to  propagate  the 
waveforms  outside  of  the  source  region  and  (2)  by  direct  calculation  using  the  same  finite-difference  code.  The 
comparison  is  very  good,  which  shows  that  the  representation  theorem  technique  is  working  correctly  and  that  we 
can  effectively  use  the  representation  theorem  to  extend  the  size  of  the  finite-difference  calculation  to  any  desired 
range,  provided  that  the  earth  can  be  modeled  as  a  plane-layered  structure  at  points  outside  the  monitoring  surface. 
The  technique  also  requires  that  the  material  properties  at  the  boundary  of  the  calculation  match  the  plane-layered 
structure.  We  have  implemented  the  technique  for  full  waveform  seismograms  and  plan  to  do  the  same  for  modal 
seismograms  and  far-field  body  waves  as  part  of  this  project.  The  reason  for  implementing  the  other  Green’s 
functions  is  that  the  modal  decomposition  allows  us  to  look  at  the  effects  on  individual  modes,  and  the  far-field  body 
wave  representation  provides  the  clearest  picture  of  the  P  and  S  waves  generated  by  the  source. 

X  1 0’’®  TRES  vs  TRES+Prose 


Figure  4.  Comparison  between  solutions  using  the  finite-difference  method  (solid  lines)  and  the  solutions 
using  the  representation  theorem  (dashed  line)  at  a  distance  of  2561  m.  The  three  curves  from 
bottom  to  top  are  the  radial,  tangential,  and  vertical  components. 

The  principal  advantage  of  this  approach  is  that  it  allows  us  to  perform  detailed  calculations  of  the  source  region  and 
then  propagate  the  results  to  distances  that  would  be  impractical  or  impossible  to  include  in  the  same  numerical 
calculation.  In  addition  to  reducing  cost  and  time,  the  hybrid  method  is  also  more  accurate,  as  numerical  dispersion 
increases  with  the  size  and  duration  of  numerical  calculations.  A  similar  approach  (in  reverse)  was  used  by  Moczo 
et  al.  (1997)  to  calculate  the  effect  of  topography  and  near-surface  structure  on  ground  motion  caused  by  a  distant 
source.  Wu  et  al.  (2004,  and  earlier  papers  starting  with  Xie  and  Lay,  1994)  used  a  similar  concept  by  performing 
near-source  finite-difference  calculations  followed  by  a  calculation  of  slowness  at  large  distances,  which  they  then 
used  to  estimate  the  energy  trapped  in  the  crust  versus  lost  to  the  mantle  for  different  scattering  scenarios.  A 
significant  advantage  of  the  representation  theorem  approach  described  here  is  that  we  can  also  calculate  waveforms 
at  any  location. 

Development  of  a  3D  Nonlinear  Finite  Element  Code 

A  difficult  task  in  this  new  project  is  development  of  a  3D  nonlinear  code  with  capabilities  similar  to  our  2D 
axisymmetric  finite-difference  code  CRAM.  CRAM  has  a  long  history  of  successful  calculation  of  explosions,  and 
has  well-tested  material  models.  It  also  has  the  capability  of  including  gravity  and  can  model  spall  and  other  free- 
surface  nonlinear  interactions.  For  the  current  project,  we  are  using  a  finite  element  code,  as  this  is  more  easily 
adapted  to  complex  3D  geometries.  In  particular,  the  cavity  can  be  any  shape,  and  the  media  can  be  nonlinear  and 
heterogeneous. 
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We  have  implemented  an  explicit  3D  Lagrangian  finite-element  algorithm  that  is  capable  of  using  multiple 
processors.  All  of  the  nonlinear  material  models  from  2D  CRAM  have  been  implemented.  The  cavity  is  placed  near 
the  center  of  the  grid  and  is  enclosed  by  a  spider  grid  that  facilitates  applying  the  pressure  boundary  condition  and 
rezoning  elements,  as  implemented  in  the  2D  axisymmetric  code,  CRAM.  Examples  of  the  spider  grid  and  outer  grid 
are  shown  in  Figure  5. 
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Figure  5.  Left:  Schematic  of  one  quadrant  of  the  spider  grid.  The  red  grid  lines  represent  the  cavity 
boundary,  the  blue  lines  represent  the  tangential  grid  lines  inside  the  spider  grid,  and  the  green 
lines  represent  the  interface  between  the  spider  grid  and  the  outer  Cartesian  grid.  Right:  The  outer 
Cartesian  grid  enclosing  the  spider  grid  (left). 

In  the  numerical  scheme,  the  wavefields  in  every  element  are  computed,  then  the  wavefields  at  each  grid  are 
computed  by  summing  the  contributions  from  each  adjacent  element  (Ma  and  Liu,  2006).  Since  the  one-point 
integration  method  is  used  for  efficiency  and  flexibility,  the  resulting  hourglass  modes  are  damped,  using  both  a 
viscous  hourglass  control  scheme  and  a  stiffness  hourglass  control  scheme  for  large  deformation.  For  grid  stability 
due  to  strong  shock-wave  propagation,  artificial  viscosity  is  also  necessary.  Both  linear  and  quadratic  viscosity 
(extended  from  2D  CRAM)  have  been  implemented. 

Benchmark  Tests 


As  a  first  benchmark  test  on  the  implementation  of  the  finite  element  method,  we  calculate  the  dynamic  response  of 
a  half-space  to  a  pressure  pulse  applied  to  the  wall  of  a  spherical  cavity  with  a  2-m  radius  buried  12  m  deep  in  a 
linear  elastic  medium  and  compare  this  with  the  analytical  results  obtained  using  the  point  source  theory  of  de  Hoop 
(Jiang  et  al.,  1994).  The  agreement  is  excellent  and  is  the  first  step  in  validating  the  implementation  of  the 
finite-element  algorithm  (Figure  6). 

As  a  second  benchmark  test,  we  calculate  the  dynamic  response  of  a  half-space  nonlinear  granite  model  (Stevens 
et  al.,  2003)  to  an  explosion  of  20  tons  yield  in  a  5-m  spherical  cavity,  buried  at  a  depth  of  102  m,  using  the  new  3D 
CRAM  and  the  existing  2D  CRAM.  Figure  7  shows  a  comparison  of  seismograms  at  the  free  surface.  The  results 
show  that  the  3D  CRAM  calculations  are  also  in  excellent  agreement  with  the  existing  2D  CRAM  calculations, 
demonstrating  that  both  codes,  though  implemented  by  different  numerical  methods,  are  capable  of  accurately 
computing  nonlinear  wave  propagation  in  rocks  due  to  explosions  in  a  cavity. 
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Figure  6.  Dynamic  response  of  a  half-space  to  a  pressure  pulse  applied  to  the  wall  of  a  buried  spherical  cavity. 
The  figure  shows  the  vertical  component  of  displacement  versus  time.  The  results  are  consistent 
with  the  analytical  solutions  (Jiang  et  al.,  1994). 
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Figure  7.  Comparisons  of  seismograms  at  three  surface  locations — 120  m  (top  left),  228  m  (top  right),  and 

330  m  (bottom).  The  top  row  shows  the  vertical  component;  the  bottom  shows  radial.  Blue  lines  are 
calculated  with  2D  CRAM;  red  lines  are  calculated  with  3D  CRAM. 
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Explosions  in  Rectangular  and  Spherical  Cavities 


To  illustrate  the  effect  of  a  nonspherical  cavity  in  three  dimensions,  we  compare  the  waveforms  from  detonating  the 
same  yield  of  20  tons  in  a  spherical  cavity  (red  in  Figure  8,  left)  and  a  cubic  cavity  (blue),  both  cavities  having  the 
same  volume.  The  recordings  at  the  4  surface  stations  are  compared  in  Figure  8  (right).  Stations  1  and  3  are  located 
at  a  distance  of  about  8-10  times  the  cavity  size  (about  3  times  the  nonlinear  region  size).  Although  the  cavity 
volume  is  the  same,  the  differences  between  the  seismograms  generated  by  the  two  cavity  shapes  are  clearly  visible 
on  the  three  components.  The  explosion  in  the  rectangular  cavity  is  found  to  cause  more  nonlinear  deformation 
around  the  cavity  and  to  generate  shear  waves  shown  in  the  tangential  direction. 
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Figure  8.  Map  view  of  cavity  and  stations  (left).  Comparison  of  velocity  seismograms  at  4  locations  (right). 
Red  and  blue  lines  correspond  to  spherical  and  rectangular  cavities,  respectively. 


CONCLUSIONS  AND  RECOMMENDATIONS 

We  are  in  the  early  stages  of  a  new  project  to  understand  3D  effects  on  seismic  radiation  from  underground  nuclear 
explosions.  We  have  partially  completed  development  of  a  3D  version  of  CRAM,  the  Lagrangian  code  we  have  used 
previously  for  performing  axisymmetric  calculations  of  underground  explosions.  We  have  also  implemented  a 
procedure  for  propagating  the  results  of  3D  source  region  calculations  to  regional  and  teleseismic  distances  using  the 
representation  theorem  to  couple  source  region  calculations  to  a  wavenumber  integration  code.  Our  plans  are  to 
complete  implementation  of  the  numerical  methods  and  then  perform  3D  calculations  to  understand  and  model  the 
effects  of  3D  source  region  heterogeneity  and  the  seismic  response  to  that  heterogeneity. 
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